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ABSTRACT 


An extension of the reacting Hz-air computer code SPARK is presented, which enables 
the code to be used on any reacting flow problem. Routines are developed that calculate 
in a general fashion, the reaction rates, and chemical Jacobians of any reacting system. In 
addition, an equilibrium routine is added so that the code will have frozen, finite rate, and 
equilibrium capabilities. The reaction rates for the species are determined from the law of 
mass action using Arrhenius expressions for the rate constants. The Jacobian routines are 
determined by numerically or analytically differentiating the law of mass action for each 
species. The equilibrium routine is based on a Gibbs free energy minimization routine. The 
routines are written in Fortran 77, with special consideration given to vectorization. Run 
times for the generalized routines are generally 20% slower than reaction specific routines. 
The numerical efficiency of the generalized analytical Jacobian, however, is nearly 300% 
better than the reaction specific numerical Jacobian used in SPARK. 

NOMENCLATURE 

Aj : reaction rate constant for jth reaction 

bj : body force of species i 

C, : concentration of species i 

Cj : time rate of change of Cj 

C p : specific heat at constant pressure 

E : total internal energy; activation energy 

‘Presently, Research Scientist working under contract for NASA Langley Research Center 


/, : mass fraction of species i 

Qj : Gibbs free energy of species i 

G r : Gibbs energy of reaction 

hi : enthalpy of species i 

h°j : base enthalpy of species i 

Kf : forward rate constant 

Kb : backward rate constant 

K eq : equilibrium constant 

Mi : molecular weight of species i 

rii : moles of species i 

nl : number of species constraint equations 

ns : number of chemical species 

nr : number of chemical reactions 

p : pressure 

q : heat flux 

R° : universal gas constant 

S : entropy 

T : temperature 

Tr : reference temperature = 298 °K 

t : time 

u : streamwise velocity 

v : transverse velocity 

Ui : streamwise diffusion velocity of species i 

i>i : transverse diffusion velocity of species i 

Vi : diffusion velocity vector of species i 

x : streamwise coordinate 

y : transverse coordinate 

Xi : mole fraction of species i 

A : second viscosity coefficient; Lagrangian Multiplier 

7 : stoichiometric coefficient; specific heat ratio 

p : density 

<7 : normal stress 

r : shear stress 

u)i : species production rate of species i 

Pi : chemical potential of species i 

p : laminar viscosity 

£ : computational streamwise coordinate 

r\ : computational transverse coordinate 

INTRODUCTION 


The SPARK [1] code is a two-dimensional Navier-Stokes and chemistry code which 

has gained wide acceptance as a means of resolving flows involving reacting hydrogen and 
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air. In addition to solving the conservation equations of mass, momentum and energy, 
it tracks nine chemical species, eighteen reaction paths, and uses finite rate chemical ki- 
netics to simulate the chemical reactions. Numerically, it is an explicit, time accurate, 
spatially elliptic finite difference code, which uses MacCormack’s technique[2] to integrate 
the equations spatially and temporally. The codes’ architecture incorporates vectorization 
whenever possible for the numerous spatial arrays in the formulation, making it highly 
efficient. One disadvantage with the present formulation of the SPARK code is that it 
is written specifically for the if 2 -a ir reaction. It would be advantageous if the reaction 
specific chemistry could easily be modified so that any chemically reacting flow could be 
modeled without substantially rewriting the code. 

For certain flows, it is informative to compare finite rate results with those obtained 
assuming frozen and equilibrium chemical approximations. The frozen and equilibrium 
limits represent the lower and upper bounds achievable from chemical interactions and 
bracket the finite rate results between them. For .ff 2 -air systems, seldom is the equilibrium 
approximation strictly valid, owing to the rapid hydrodynamic timescales in the problems. 
It is, however, useful to perform bracket analyses on reacting flows, to see how nearly frozen 
or equilibrated each flow is. Nozzles have traditionally been calculated assuming the frozen 
and equilibrated chemistry approximations, thus relying heavily on the equilibrated flow 
capabilities. 

For flows where the chemistry timescale is orders of magnitude faster than the hydrody- 
namic timescales, the equilibrium approximation becomes a realistic means of accounting 
for the chemical relaxation processes, at considerable savings over full finite rate calculi 
tions. In particular, high temperature air chemistry fits into this category, and could be 
reasonably approximated with an equilibrium chemical assumption. 

The most efficient means of calculating a particular flow depends in the ratios between 
the chemical and hydrodynamic timescales. For fast (slow) chemistry, the equilibrium 
(frozen) approximation can be used effectively. If the timescales are of the same order, a 
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finite rate calculation must be used. Which chemical species and reaction paths are needed 
to resolve the physics, or whether a frozen or equilibrated flow is applicable, is generally 
not known before starting the calculation, and must be iterated on to find the most efficient 
method for each problem. This iterative process makes a generalized SPARK code very 
desirable. 

The principle motivation for this work is the need for a modified version of the existing 
SPARK [l] code, capable of accommodating any reaction system. It must be easy and 
systematic to modify, thus reducing programming time and errors. The SPARK code 
is presently capable of tracking the governing fluid mechanical equations, as well as any 
number of chemical species. Since the chemistry portion of the SPARK code is modular 
in structure, the underlying framework need not be changed to accommodate different 
chemical species, reactions or reaction rates. Only the subroutines in which these features 
are calculated must be modified. 

Specifically, the SPARK code is modified in the following way. A generalized properties 
routine is added to the code, in which the thermodynamic and kinematic properties of 
the reacting species are introduced. Where or how these properties are obtained from 
the literature is left up to the user. The existing reaction rate routine is replaced by a 
generalized routine, capable of determining the reaction rate of any reacting system. The 
chemical Jacobian routine which calculates the numerical Jacobian of the H 2 -& ir reacting 
system is replaced by a generalized numerical Jacobian routine. In addition, a generalized 
analytical chemical Jacobian is developed to promote numerical accuracy and efficiency. 
Finally, a routine which calculates the equilibrium species concentrations given a particular 
temperature and pressure is added. All subroutines which are developed are built around 
the existing SPARK code structure in which the unknowns (for instance the species reaction 
rates) are calculated simultaneously at all grid points in the domain. As in SPARK, a 
judicious attempt is made to utilize vectorization whenever possible. 

The usefulness and need for a generalized reaction scheme as well as a generalized Ja- 
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cobian for the SPARK code is clearly evident from the above discussion. This code could 
be used during an optimization study to vary reacting species and paths systematically, 
or to model any chemically reacting system. Thus, given sufficient knowledge of the rel- 
evant reaction species, the appropriate reaction routes, and accurate reaction rates, any 
chemically reaction flow can be resolved in a spatially and temporally developing manner. 
Although the code was designed specifically for the SPARK code, it is hoped that the basic 
ideas can be extended to other chemistry codes which need similar modifications. 

ANALYSIS 


Governing Equations 

Since extensive detail of the SPARK code exists in Reference [1], particularly the 
governing equations and numerical techniques needed to solve these equation, only a brief 
summary of the equations will be presented here. The chemistry equations will be ad- 
dressed with more rigor since they are the emphasis of this work. 

The equations governing the conservation of mass, momentum in two dimensions, and 
energy, and the species conservation equations can be written as [1] 

Pt + (pu) x + ( pv) y = 0 (1) 

ns 

C pu)t + ( puu - a x ) x + (puv - T xy ) y = pJ2 fib* ( 2 ) 

1=1 

ns 

(pv)t + ( puv - Ty X ) x + (pvv - o y ) y = pJ2 fihiy ( 3 ) 

t=l 

ns 

(pE) t + \{pE - o x )u - t xv v + q x \ x + \(pE - (T v )v - r yx u + q y ] y = P^fM • {V + V) (4) 

1 = 1 

(pfi)t + \p{u + u)fi] x + [p(v + v)fi] y = u (5) 

Ox = -p + + v v) + 2 A tu ® 
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( 6 ) 


(7) 


<7y = — P + A (Us + Vy) + 2flVy 


Qx 


T X y = 7y X = + V s ) 


ns ns ns V.1T) 

= —kT x + P £ hifiUi + i?°T Z) ~ «j) 


t=l 


^ ^ MD 
•= 1 i=i iVl ' u ') 


ns ns ns y r\ 

«, = -*r, + />E *■ /. ». + *°r X) £ inf ft - s, ) 

i=l »=1 , = 1 W1 i U H 


F sr hf p.J + v 2 


( 8 ) 

(9) 

( 10 ) 

( 11 ) 


h, = /t,° + f T C p .dT (12) 

JT r 

% = A,- + B,T + C.T 2 + AT 3 + AT 4 (13) 

R 

ns r 

P = pR°T'£-£- 

U M ‘ 

Equations 1-4 represent the continuity, u- and v- momentum, and energy equations. 
Equation 5 is the individual species mass conservation equation. Equations 6-13 describe 
how particular terms in equations 1-5 are modeled. Note that the enthalpy of the mixture 
is the summation of the partial enthalpies of the chemical species, and that C p is modeled 
with a 4th order polynomial [3]. 

Generalized Chemistry 

Chemical reaction rate expressions are usually determined by summing the contribu- 
tions from each relevant reaction route (or path), to obtain the total rate of change of 
each species. Each path is governed by a law of mass action expression in which the rate 
constants can be determined from a temperature dependent Arrhenius expression. The 
reaction rate terms u> present in equation 5, written in this form can be expressed as 
follows: 
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3 = 1 ,nr 


(15) 


ns Kf . ns 

£%•<<?, ^ £ 7 ;,c, 

»=1 t=l 


ns f 


u i = M i C t = M i '£.('1i<--ij(Wl, n c ” - •'f*, II «“1 i = l,ns (16) 

j—\ m= 1 m=l 

where equation 15 describes the relevant chemical reaction paths for the system of species 
present in equation 5, and equation 16 is the reaction rates for each species i, as determined 
by the law of mass action. 

The reaction constants Kj i and Kb } . are determined as follows: 

K U = A,T"t exp(-j^) (17) 


X. = 




(18) 


K = f 1 l An cxpf ^ Gr > ) 

A *9y l Ro rp) Pi R o T ) 


(19) 


A ° R ; = 2 J = 1, nr 


1=1 


1=1 


( 20 ) 


§ = - r in t) + Ar* + §r s + + ^r 1 + f,- - c,r (21) 

The forward reaction rate for each reaction path is determined from equation 17 from 
an Arrhenius law expression. The reverse reaction rate is back-calculated in equation 18 
from the equilibrium constant defined in equation 19. Equations 20 and 21 describe how 
the change in Gibbs free energy for each reaction is calculated. 

Generalized Jacobian 

For some reacting flows, the chemical Jacobians of the chemical system must be evalu- 
ated. Certain chemical systems are numerically unstable when integrated with a time step 

comparable with the CFL condition. One remedy for this instability is treating implicitly 
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the destabilizing chemical source terms u>i [10]. This technique requires the Newton lin- 
earization of the chemical source terms about the previous time step and results in the need 
for a chemical Jacobian relating the rate of change of each source term with respect to each 
species. This Jacobian can be calculate either numerically or analytically, and presents no 
real problem once a relevant system of reactions and species has been introduced. 

Analytically, the chemical Jacobian is obtained by differentiating equation 16, and is 
defined by: 


dCi 

dC m 


nr TS ns i fr ns » 

= E(%-< - v<)[-P= n n <?'■'] •• = i,»» 


( 22 ) 


i=i i=i i=i 

The numerical Jacobian is calculated by approximating the limit formula for the deriva- 
tive as: 


dCi C i {P,T,C u ..C m - 1 ,C m + AC m ,C m+u ..C ns ) 

dC m a&o A C m [ } 

Each species concentration is individually perturbed while holding all others constant. 

The change in reaction rate is then assessed. Note that, in both Jacobian routines, the 

third body concentrations relevant in some of the chemical reactions are assumed constant 

to simplify the procedures. This is a valid assumption since the perturbation concentration 

is negligible compared with the combined concentrations of all species present. Also note 

that the third body collision efficiency is assumed equal to one for all species. This has 

been shown to be an invalid approximation [9] under certain flow situations, but is used 

here for simplicity. 

Equilibrium Equations 

Calculations which assume the flow to be in thermodynamic equilibrium, require an 
additional equation and do not rely on equation 5 for the local species mass fractions. 
The equilibrium state is characterized by the fact that the total Gibbs free energy of the 

mixture is at a minimum at this point, subject to the constraints of conservation of atomic 
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species [5]: 


ns 

G = YZ 

(24) 

i=i 


ns 


Y2 a ij n j — b° = 0 / = 1, nl 

(25) 


;'=i 


Equation 24 defines the Gibbs free energy of the system, defined as the summation of the 
individual components free energy. Equation 25 relates the species constraint equations 
dictating the conservation of atomic mass of each elemental species. These equations 
are most easily solved using the method of Lagrangian multipliers. With this technique 
equations 24-25 become: 

ns 

(Aj + ) ' XfCL^j = 0 

«=1 

nt 

YY aijrij - b° = 0 / = 1 , nl (26) 

i=i 

where the A,- are the Lagrangian multipliers, and nl is the number of constraints. Note 
that equation 26 is not solved when finite rate chemistry is being used. It is introduced 
only for the purpose of enabling SPARK to be run in the equilibrium mode. 

GENERALIZED REACTION PROGRAM 

The generalized chemistry routine consists of five subroutines: 1) PROPS, 2) GEN- 
RATE, 3) GENJACA, 4) GENJAC, and 5) EQUIL. Subroutine PROPS specifies the con- 
stants for the chemical species begin studied, so that its exact form depends on the chemical 
system being modeled. The other four subroutines have been written in a general form 
which applies for any chemical system obeying equations 1 - 21. 

In subroutine PROPS, all of the constants used to define the chemical species’ thermo- 
dynamic and kinematic properties, as well as the appropriate reaction rate constants are 
specified. The exact structure of PROPS is strongly dependent on the chemical system 

being simulated, on how the thermodynamic and kinematic features are modeled, and 
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what form of chemical kinetics is used to calculate reaction rates. Implicit, however, is the 
assumption that all chemical systems behave according to the law of mass action through 
equations 15-23. 

For SPARK, nine chemical species are used to represent the H<i - air chemical reaction 
system. The thermodynamic properties (enthalpies, entropies, Gibbs free energies, and 
specific heats) are specified in PROPS from fourth order polynomial fits obtained from 
McBride [3]. The species’ kinematic properties are specified from analytic expressions 
obtained from Williams [4] and White [11]. Specifically, the molecular viscosities and 
conductivities are modeled using Sutherlands law [3,11,13], and the binary diffusivities are 
obtained from Williams, and White [4,11]. Arrhenius type expressions are used for the 
reaction rates (equations 15 - 21), and appropriate reactions and coefficients are obtained 
from Jachimowski, and Dixon-Lewis [14,15]. The exact form of the properties routine will 
certainly change with different chemical reactions and is only presented here as an example. 

GENRATE 

Subroutine GENRATE calculates the reaction rates u> for each of the chemical species 
as a function of space and time. The control elements of the routine are two (nr x ns+1) 
bit vectors describing and , the stoichiometric coefficients i for each reaction j. If qfj. 
is zero, then species i does not occur in reaction j. If qr- is greater than one, then species i 
occurs 7,-y times in reaction j. Note that there are ns+1 apparent species in the vector. The 
last species is the third body species, defined as the summation of all individual species 
concentrations present in the mixture. The third body becomes important in some reaction 
paths where three body collisions are required to stabilize the reaction. 

The routine is essentially composed of two parts. In the first part, all the forward 
and backward reaction rates are calculated from equations 17-19, as well as all the species 
concentrations. The second portion of the routine calculates the species reactions rates 
directly from equation 16. This is done by summing the individual contributions to the 
species creation and depletion rates over all the reaction paths. Because the GENRATE 
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routine is completely general it only depends on properties provided to it from the PROPS 
routine. 

GENJACA 

Subroutine GENJACA calculates the Jacobian matrix J = analytically. Like 

the GENRATE routine, its control elements are the two bit vectors and 7L describing 

which species i are involved in reaction j. The structure of the code is nearly identical to 
the GENRATE routine except that equation 22 is used instead of equation 16. The first 
part of the code calculates all the forward and backward reaction rates from equations 17 
- 19 and the species concentrations. The second portion of the routine uses equation 22 
to calculate each component of the (ns-1 x ns-1) Jacobian matrix describing the rate of 
change of reaction rate i with respect to species m. Again, as in the GENRATE code, 
each element of the Jacobian is the summation over all the reactions of the forward and 
backward contribution to the derivative. Also note that no where does division by species 
” C m ” occur in the program; rather the exponent 7 i} - is decreased by one to account for this 
effect. The GENJAC program is completely general, and depends only on the property 
information given to it by subroutine PROPS. 

GENJAC 

Subroutine GENJAC calculates the Jacobian matrix J = numerically. It approx- 
imates each element of the matrix according to equation 23 by perturbing each species 
independently, and accessing the change in the reaction rate. Care must be take to en- 
sure that the perturbation is of the proper order. If the perturbation concentration is too 
small, its effect may be indistinguishable from machine roundoff error. If it is too large, 
the approximation to the derivative may be poor. To avoid these problems, a median 
perturbation is used in which the perturbed value of the concentration is approximately 
equal to the change experienced by that concentration over the previous time step, or 10 -6 , 
which ever is larger. 

The routine is again composed of two parts. In the first, all the reaction rates are 


11 


determined from equations 17 - 19 as well as the unperturbed species concentrations. 
Then a perturbation concentration is proposed and the reaction rates for each species are 
determined in a manner identical to that used in GENRATE: i.e. the cumulative addition 
to each species forward and backward reaction rate from each reaction route. Finally, the 
simple calculation of [u),(C m + AC m ) -Ui(C m )]/AC m yields the explicit values for J. Again 
the routine is completely general and only depends on the properties given to it by the 
PROPS routine. 

EQUIL 

The EQUIL routine calculates the species concentration profiles at each spatial location, 
assuming the species to be in chemical equilibrium. The routine is based on a scalar version 
of an equilibrium routine developed by Gordon and McBride [5,6]. Since details about 
the theory and implementation of the Gordon and McBride codes can be found in the 
literature, only a brief review of the theory will be presented here. 

The EQUIL routine solves equation 24 at each gridpoint. Since the equation is a 
nonlinear algebraic system, a Newton iteration is used to seek the root. The iteration is 
relaxed initially to ensure more robust convergence properties. The relaxation is deter- 
mined by monitoring the rate of change in iteration space of the species concentrations. As 
the solution approaches the root, the relaxation is eliminated, thus resulting in quadratic 
convergence to the root. Generally, between 5 and 20 iterations are required to drive the 
residual down 10 orders of magnitude (far more than is actually needed), depending on 
the accuracy of the initial guess. 

The routine minimizes the total Gibbs free energy of the mixture using the method of 
Lagrangian multipliers. The species Gibbs free energy is calculated from the thermody- 
namic data available in the PROPS routine using equations 21. The constraints on the 
atomic mass b° are imposed by summing the initial levels of atomic species present at the 
spatial location. No provision is made for spatially varying constraint equations, since 
the computational efficiency of this technique was found to be uncompetitive with simple 
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finite rate calculations. The minimization assumes that the pressure and the enthalpy of 
the mixture remain constant during the minimization. Thus, new species concentrations 
and a new temperature are the resultant of the minimization. 

The coupling between the SPARK code and the EQUIL routine is very arbitrary. For 
most time steps, the SPARK code is run in frozen mode, and the frozen species values 
are updated periodically. This eliminates the need to explicitly account for the species 
continuity equations, and makes the procedure only slightly less efficient than a frozen 
flow calculation. However, the entire calculation is no longer time accurate due to the 
lagging in the chemistry. The species concentrations are assumed frozen between each 
call. The call updates the local species concentrations and temperature subject to the 
constraints of mass conservation, and constant pressure and mixture enthalpy. This was 
found by Davy [8] to be the most efficient method of coupling the hydrodynamic and 
species equations. 

For each problem the most efficient method of updating can be found by trial and error. 
In general, the EQUIL routine should be called as infrequently as possible because of the 
amount of CPU time required by the Newton iteration. The problem, geometry, initial 
conditions, and various other parameters dictate this coupling. In general, the EQUIL 
routine is called about every 25-50 iterations. 

The EQUIL routine vectorizes in the spatial direction, thus promoting computational 
efficiency. In addition, the EQUIL routine is written to accommodate the need (the SPARK 
code is presently not formulated with this feature) for spatially varying constraint equa- 
tions, such as would be the case in diffusionally driven reacting flows. For these problems, 
the situation exists where certain portions of the flow will converge more rapidly than 
others. In this case, the converged elements of the vector are eliminated and the vector is 
repacked to again ensure vector efficiency. 

RESULTS 

The generalized reaction routines have been run in conjunction with the SPARK code 
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to access their validity, accuracy and efficiency. The code exhibits striking simplicity in 
changing from one chemical reaction system to another. If the thermodynamic and kine- 
matic information is readily available for the reacting species, testing various reaction paths 
is a matter of forming the control bit vectors 7^ and 7?., and increasing the accompanying 
vector lengths. An obvious advantage in the formulation is that the chances for human 
error are greatly reduced. For example, changing the code from a Carbon based reacting 
flow to a Hi-dAr reacting system can be done in a matter of minutes if the properties of all 
the relevant species are presently known. 

The scheme exhibits two minor disadvantages; 1) increased program complexity, and 
2) decreased computational efficiency. A direct comparison between the existing SPARK 
reaction specific rate routine and the generalized scheme, each involving nine species and 18 
reactions, shows the generalized scheme to require 20% more CPU time. A 30% increase 
in CPU time is required when using the generalized numerical Jacobian to replace the 
reaction specific Jacobian in SPARK. Because only a portion of the total computer time 
used by a code is used in these subroutines, the total increase in run time for SPARK 
is a somewhat smaller percentage. For example, for typical fully explicit runs, the 20% 
increase in the generalized reaction routine translates to about a 10% increase in total run 
time (50% of total time is spent in the reaction routine for a time accurate calculation). 
For implicit calculations where both the generalized reaction and Jacobian routines are 
used, total computer time increases greater than 20% are found since nearly 65% of all the 
time is spent in calculating the numerical Jacobians at each time step. 

In comparing the run times of the generalized analytic and numerical Jacobians (a reac- 
tion specific analytical Jacobian does not exist at present in the spark code), the numerical 
Jacobian is about 285% more costly. Comparing the generalized analytic Jacobian with the 
reaction specific numerical Jacobian finds the reaction specific Jacobian over 200% more 
costly. Thus, the use of the generalized analytic Jacobian in SPARK decreases run-times 
by 25% for fully implicit runs! A reaction specific analytical Jacobian could be expected 
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to increase the efficiency of fully implicit SPARK by approximately 30% based on these 
findings. 

The use of the EQUIL routine increases the run-time of the SPARK code about 10% 
over the frozen run-time. It should be noted that the frozen and the equilibrium chemistry 
modes for the SPARK code are at least an order of magnitude more efficient than finite 
rate results using full chemistry. 

Figure 1 shows the final calculations performed during the testing of the equilibrium 
routine. In each figure the frozen and equilibrium limits were found for a chemically 
reacting system. Then the reacting rates were increased and decreased from their actual 
values to see if in the limits the frozen and equilibrium values were approached. This 
was achieved by multiplying the pre-exponential factors for all the reaction paths by an 
arbitrary constant. By varying this constant, the chemistry time could be varied from 
” infinitely fast” to ” infinitely slow”. Figure 1 shows this test on the H 2 -&. 1 T system. The 
thrust levels from the nozzle are plotted against the log of the arbitrary multiplicative 
constant. The figure shows that asymptotic convergence is achieved at both the fast and 
slow extremes. 

Figures 2-11 are presented to exhibit the capabilities of the generalized version of 
the SPARK code. They are not meant to display intricate details of specific reacting 
flows, but rather show a situation where a generalized version of the SPARK code is very 
convenient. The figures show the flow through a supersonic nozzle as calculated assuming 
two different chemical reaction systems. The question being addressed here is the overall 
importance of reacting nitrogen in nozzles of similar design and inlet conditions. This flow 
is complicated by the stiff nature of the reaction terms near the inlet [7] making the use of 
a chemical Jacobian essential. The nozzle inlet conditions are representative of combustor 
exit conditions in which H 2 and air were combusted at stoichiometric conditions, and are 
very nearly in chemical equilibrium. The two reacting systems, being compared here are, 
1) 9 species and 18 reactions assuming nitrogen to be inert, and 2) 15 species and 35 
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reactions accounting for nitrogen chemistry. 

Figures 2-5 show the velocity vectors, the pressure contours, the temperature contours, 
and the water mass fraction contours, respectively, for a two-dimensional supersonic nozzle. 
The lower surface represents the lower wall of the nozzle, and the upper surface is the 
centerline, or line of symmetry of the nozzle. The nozzle is calculated on a 30 by 20 grid, 
which by no means fully resolves the boundary layer along the wall. Figures 6-11 show a 
comparison of temperature profiles as calculated assuming various chemistry models. All 
of the figures show the axial variation of the temperature along the wall, J=l, midpoint, 
J=ll, and the centerline, J=20. Figures 6-8 show a frozen, a finite rate chemistry, and 
an equilibrium calculation, calculated assuming nitrogen to be inert, with the nine species 
reaction scheme. Figures 9-11 show the flows calculated in the same manner with full 
account of any nitrogen reactions which may be present. Comparisons of the differences 
resulting from each chemistry model are easily made using the generalized code. Thus, 
an appropriate determination of the best chemistry model to be used for this flow can be 
made. 

Presently, work is continuing to optimize the generalized routines’ numerics and to 
make the properties more applicable to non-SPARK type CFD codes. 

CONCLUDING REMARKS 

A generalized reaction code has been developed as an add-on feature for the presently 
existing SPARK code. It has been written so as to eliminate the tedious steps necessary in 
converting the SPARK code from one particular reacting system to another. Its usefulness 
has been demonstrated in modeling a variety of air chemical reaction systems. The 
generalized routines have a slight decrease in efficiency associated with them, so that if 
numerous runs are to be made with one reaction system it is advantageous to develop a 
reaction specific code. The study also dictates that the use of an analytic Jacobian should 
be made whenever possible. 
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H2 - AIR BRACKET STUDY 
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Figure 1: air asymptotic finite rate chemistry study. Chemistry approaches the frozen 

and equilibrium limits as the rate constants are adjusted. 



Figure 2: Velocity vectors for viscous flow through nozzle. 



Figure 3: Temperature contours for viscous flow through nozzle. 


Figure 4: Pressure contours for viscous flow through nozzle. 
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Figure 5: Water mass fraction contours for viscous flow through nozzle. 
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Figure 6: Frozen temperatures along the wall, midpoint and centerline. 
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Figure 7: Reacting temperatures along the wall, midpoint and centerline. 
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Figure 8: Equilibrium temperatures along the wall, midpoint and centerline. 
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Figure 9: NO x Frozen temperatures along the wall, midpoint and centerline. 
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Figure 10: NO x Reacting temperatures along the wall, midpoint and centerline. 
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Figure 11: NO x Equilibrium temperatures along the wall, midpoint and centerline. 
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